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We report results of the study of solitons in a system of two nonlinear- Schrodinger (NLS) equations 
coupled by the XPM interaction, which models the co-propagation of two waves in metamaterials 
(MMs). The same model applies to photonic crystals (PCs), as well as to ordinary optical fibers, close 
to the zero-dispersion point. A peculiarity of the system is a small positive or negative value of the 
relative group-velocity dispersion (GVD) coefficient in one equation, assuming that the dispersion is 
anomalous in the other. In contrast to earlier studied systems of nonlinearly coupled NLS equations 
with equal GVD coefficients, which generate only simple single-peak solitons, the present model gives 
rise to families of solitons with complex shapes, which feature extended oscillatory tails and/or a 
double-peak structure at the center. Regions of existence are identified for single- and double-peak 
bimodal solitons, demonstrating a broad bistability in the system. Behind the existence border, 

Othey degenerate into single-component solutions. Direct simulations demonstrate stability of the 
solitons in the entire existence regions. Effects of the group- velocity mismatch (GVM) and optical 
loss are considered too. It is demonstrated that the solitons can be stabilized against the GVM by 
means of the respective "management" scheme. Under the action of the loss, complex shapes of the 
solitons degenerate into simple ones, but periodic compensation of the loss supports the complexity. 

Q_| ' PACS numbers: 42.70.Mp; 42.70.Qs; 42.81. Dp; 05.45.Yv 

I. INTRODUCTION AND THE MODEL 

Theoretical and experimental studies of various artificial optical media have recently drawn a great deal of interest, 
see, e.g., Refs. [l[ and references therein. Among these, a well-known class is formed by metamaterials (MMs) 
based on the intrinsic periodic structures, built with subwavclcngth characteristic scales, that feature a combination 
\& , of negative dielectric permittivity and magnetic permeability, and thus give rise to the negative refractive index 2] 
MMs promise a number of potential applications impossible in ordinary optical media, such as superlensing d _ 
MMs are also described as "left-handed" waveguides, as they were originally predicted as those where the wave vector 
of the electromagnetic wave is antiparallel to the usual right-handed cross product of the electric and magnetic fields 
d . The intensive work has resulted in the creation of MMs featuring the negative refraction at optical frequencies [H . 
Existing MMs are usually assembled as periodic arrays of metallic split-ring resonators, with various modifications of 
this basic setting. 

Theoretical analysis of various nonlinear effects in models of MMs including cubic 

d~d or quadratic [13, El 

terms 

has also drawn considerable attention. Some of these studies predict the existence of solitons in models of the MM type 
In particular, a theoretical model for the XPM (cross-phase- modulation)-mediated interaction of co- propagating 
waves, carried through the MM at different frequencies, is similar to the well-known two-wave model of ordinary 
media (where, in particular, it may give rise to the modulational instability in the case of the normal group-velocity 
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dispersion (GVD) |13|, and to domain- wall patterns [14(, in the same case). However, an essential peculiarity of the 
MM model is that the co-propagating waves may have widely different GVD coefficients, even with opposite signs, 
corresponding to normal and anomalous dispersion d- Another type of artificially built optical media where this 
feature may be realized is represented by photonic crystals (PCs) and PC fibers [l5|, [l6|]. In fact, these possibilities 
illustrate the versatilely of artificial optical systems, which allow one to engineer required properties, that may often 
be made very unusual. 

The main subject of the present work is constructing families of stable solitons within the framework of the system of 
XPM-coupled NLS (nonlincar-Schrodinger) equations describing the co-propagation of two waves with local amplitudes 
u(z, t) and v (z, r), carried by different frequencies: 

iu z + (l/2)u TT + (\u\ 2 +2\v\ 2 )u = 0, (1) 
iv z +icv T + (l/2)D 2 u TT + (\v\ 2 +2\u\ 2 )v = 0. (2) 

Here z and r are, as usual, the propagation distance and reduced time, the GVD coefficient in Eq. ([1]) is normalized 
to be D\ = 1 (which implies the anomalous sign of the GVD for wave it), D 2 , that may have either sign, is the relative 
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dispersion for the second wave (D 2 < implies that wave v propagates with normal GVD), c is the walk-off parameter 
(alias the group- velocity mismatch, GVM, between the two carrier frequencies), which is an essential ingredient of the 
respective MM model p}, and it is assumed that effective Kerr-nonlinearity coefficients are equal for both waves. 

It is relevant to mention that the GVM between the fundamental-frequency and second-harmonic waves also plays 
an important in role in the MM model with the \^ nonlinearity. An interesting feature of that system is that the 
phase locking of the second harmonic to the fundamental one, which is strongly affected by the GVM, may effectively 
impose the action of the negative refractive index on the second harmonic even in the situation when it directly affects 
only the fundamental frequency [Il| . 

Stationary solutions to Eqs. JTJ) and ((2|) with two independent wavenumbers of the components, K u>v , are looked 
for as 



In the case of c = 0, the stationary waveforms U, V may be assumed real, obeying the equations following from the 
substitution of expressions ((3]) into Eqs. (JJ) and ([2]): 



of strong losses is an inherent property of MMs [I7[, and, moreover, it has been rigorously proved that the negative 
refraction is not possible in a lossless medium fijj. Nevertheless, schemes have been proposed to essentially compensate 
the MM loss, using the matched impedance [U, or various gain mechanisms [20| . In any case, it makes sense to 
study solitons in conservative models of the MM type [1, Q . In addition, the model based on Eq s. p] ) and © may 
also be realized in PCs and PC fibers, where losses are present too, but as a less severe problem [I5(. On the other 
hand, as concerns solitons, a relevant approach may be to include a lossy MM sam ple , along with an amplifier, as 
elements into an optical cavity, and predict cavity solitons possible in such a setting [2lj . 




As mentioned above, Eqs. |T]) and @ also apply to the description of the co-propagation of XPM-interacting waves 
in usual optical media, such as optical fibers [I3[ , although in that case most works have been dealing with the case 
of equal GVD coefficient, D 2 = 1. Nevertheless, the case of D 2 < is relevant too, if the two carrier waves are 
chosen at wavelengths placed at opposite sides of the zero-dispersion point of the optical fiber. Moreover, a matched 
pair of the wavelengths may be chosen, so as to nullify the GVM between them (c = 0). The latter setting may 
find applications to fiber-optic telecommunications, as the normal-GVD wave, v (with D 2 < 0) may carry a stable 
periodically modulated wave, that can be used as a support structure suppressing the jitter in the data-carrying 
soliton stream launched in the anomalous-GVD mode [22j |. 

In previous studies of solitons supported by the XPM-coupled NLS equations, families of bimodal (two-component) 
soliton solutions were studied in detail, using both numerical methods (23l . |24| and the variational approximation [24| , 
but only for the case of D 2 — I and c = 0. In the general case, those solitons may feature the "elliptic polarization", 
i.e., an arbitrary ratio of energies in the two components, 



(in other words, the solitons may exist with arbitrary positive values of ratio R = K v /K u of the two wavenumbers), 
but they were simplest fundamental single-humped solitons. The coupled-NLS system with D 2 = 1 may have multi- 
hump soliton solutions too, but they are expected to be unstable [25l [26| . In this paper, we aim to demonstrate 
that, for values D 2 < 1, including small negative values of the relative GVD coefficient (i.e., the case of weak normal 
GVD in the ^-component), Eqs. {J), © and (|4]) give rise to very different families of stable single- and multi- humped 
solitons, which, in particular, entails bistability of the soliton states in a broad range of parameters. A characteristic 
feature of the solitons revealed by the analysis at small positive and negative values of D 2 is the existence of weakly 
localized oscillating tails attached to the main "body" of the bimodal soliton. 

It is relevant to mention that multi-hump solitons in multi-mode optical systems were first predicted (27j and exper- 
imentally observed [2a ] in photorefractive media, with the saturable nonlinearity, which induces the XPM interactions 
between different guided modes. The stability of such multi-hump states was later demonstrated in a rigorous form 
[H, [2t| . These studies were also extended into the two-dimensional setting [3(| • 

Basic results concerning the existence and stability (including the bistability) of the solitons in the present model 
without the GVM (c = 0) are summarized in Section 2. A the end of this section, we consider the action of the 
loss on solitons, and periodic compensation of the loss. Section 3 is dealing with an approach that allows one to 
stabilize solitons against the GVM by means of the "mismatch-management" technique: we consider a model with 



{u (z, t) , v (z, t)} = exp (iK UiV z) {U (r), V(r)} . 



(3) 





E V J±™ Mr)] 2 dr 



(5) 
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coefficient c which periodically jumps between large positive and negative values, making the average GVM equal to 
zero. Physically, this implies the use of a layered medium composed of segments with opposite values of the GVM. In 
fact, the GVM is known to be a serious issue in a different context, viz., an obstacle to the creation of temporal solitons 
in optical media with the quadratic (x^) nonlinearity [3l|, H^]. In that context, the GVM-management technique 
was recently elaborated in a theoretical form [331 ]. Actually, it resembles the previously proposed "tandem" schemes, 
which provide for an effective compensation of the GVM by way of a periodic alternation of linear and x^-nonlincar 
layers [34j. The paper is concluded by Section 4. 



II. SOLITONS AND THEIR STABILITY REGIONS 



A. The numerical method 



The underlying equations (UJ, ([2]), as well as their stationary counterparts, Eqs. ((4]), were solved in the numerical 
form by means of a pseudospectral method, which was adjusted to the present model, following the lines of work s 1351 ] 
(the application of this method to soliton solutions of equations of the NLS type was recently elaborated in Ref. [361]). 
The integration domain of variable t, of width T = 40, was sufficient to completely study the shape and dynamics 
of individual solitons, including those which feature extended "tails" (simulations of interactions between solitons 
might require using a larger domain). The domain was divided into N sub-domains, the fields in each one being 
approximated by truncated power expansions, 

M / r _ T (0)\ fc 

{uj (t,z), Vj (r, z) } = ^2 {u jk (z) , v jk (z) } ( ^ \ , (6) 

j = 1, N, where rj ^ is the midpoint of the sub-domain, and At = Tj (2N) its half-width. The substitution of 
expressions ([5]) into Eqs. UJ, @, truncation of ensuing expansions, and the use of the collocation at Chebyshev points 
and the continuity conditions for the functions and their first derivatives at boundaries between the sub-domains lead 
to a system of ordinary differential equations (ODEs) for amplitudes Ujk (z) and Vjk (z), which were numerically solved 
by means of an unconditionally stable implicit Crank-Nicolson scheme. Accordingly, stationary solitons were looked 
for as z-independent solutions of the ODE system, reducing it to a system of algebraic equations. The latter one was 
solved using the Newton's method. Sufficient accuracy of the results could be usually achieved with M = N = 8. 

As we demonstrate below, the choice of the pseudospectral method is essential for obtaining reliable numerical 
solutions, while a more common approach, based on the split-step fast-Fourier-transform (SSFFT) scheme (see, e.g., 
[3?|). may encounter problems in the most essential case considered in this work, viz., simulations of solitons with 
well-pronounced "tails" . Further details of the numerical procedure, which may be useful for other applications too, 
will be reported elsewhere. 



B. Two types of solitons 



The numerical solution produces stationary solitons of two essentially different types, viz., ones with a single peak 
(hump) at the center of the soliton, as shown in Fig. [TJ and double-hump solitons with two separated main symmetric 
peaks, see an example in Fig. [5J These figures also illustrate the bistability supported by the present model, as the 
single- and double-hump solitons displayed in Figs. QJa) and [2] exist at exactly the same value (E = 5) of the total 
energy, 



E = 



-foe 

' u(t)\ 2 + \v(t)\ 2 dr (7) 



(cf. Eq. ||SJ)), and share a common value (R = 1.13) of the wavenumber ratio that was defined above, R = K v /K u . 

Further, Fig. [3] shows that, in either case (single- or double-peak structure), the solitons develop slowly decaying 
oscillatory tails as the relative GVD coefficient in Eq. D2, takes smaller values. In fact, the problem posed by 
strong tails in bimodal systems with (slightly) normal GVD in one component is well known in the studies of optical 
models with the x^ nonlinearity, as in the cases of practical interest x^ crystals feature normal GVD at the second 
harmonic [3l|. While, strictly speaking, the tail in the normal-GVD component cannot vanish at r — * ±00, it has been 
concluded that one can identify a well-defined region in the respective parameter space, including a finite interval of 
negative values of D 2 (in the present notation), where the amplitude of the normal-GVD field assumes so small values 
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(a) (b) 

FIG. 1: Examples of stable bimodal solitons with the single central peak. The total energy of the soliton, defined as per Eq. (0, 
is E = 5 in (a) and E = 2 in (b), respectively. The other parameters are D2 = 0.37, c = 0, and the ratio of the wavenumbers 
of the two components is R = K v /K u = 1.13 in both cases. 
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FIG. 2: An example of a stable double-hump soliton found at the same values of parameters as in Fig. HI a), including the total 
energy and the wavenumber ratio. 

at r — > 00 that the nonvanishing portion of the tail is invisible and may be regarded as nonexistent, for all practical 
purposes. A transition from this case to the situation with a conspicuous long-range tail is sharp, and may be easily 
identified [32tl33l. l38j . In the present model, the tail is also, strictly speaking, nonvanishing in the case shown in Fig. 
GJ]Jb), which pertains to D2 = —0.3. However, the figure demonstrates that, in practical terms, the tail is absent at 
large values of |r|. 

Note that Fig. [3jb), which features, in addition to the well-developed tails, a slightly split peak in the v component, 
also illustrates the transition between the single-hump and dual-hump types of the solitonic shape. Particular values 
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(a) (b) 

FIG. 3: Examples of stable stationary solitons demonstrating the growth of the "tails" with the decrease of the GVD coefficient 
in the second component, D2, at c = 0, E = 2, R — 2. In (a) and (b), D2 = 0.1 and D2 = —0.3, respectively. 

of the scaled parameters for which Figs. [l][3] have been plotted (as indicated in captions to the figures) were taken 
pursuant to typical values of physical parameters in the lossless Drude model, which was used in Ref. [§| to derive 
the coupled system of Eqs. (p} and (|2|) for the bimodal propagation of electromagnetic waves in the region of the 
negative refractive index. 

Unlike the situation for simple fundamental solitons which was considered for D2 = 1 [H| j the variational approxi- 
mation (VA) for solitons is not helpful in the present case, as one would have to use a complex ansatz approximating 
the oscillatory tails, as well the possibility of having the split peak at the center of the soliton. The use of the VA 
for the description of "tailed" gap solitons in the model with the self-defocusing nonlinearity and periodic potential 
(optical lattice) [39[ , as well as for the transition from unsplit to split solitons profiles in two-component models , 
has demonstrated that this may be possible in a limited range of parameters, but the necessary calculations are quite 
cumbersome. Actually, direct numerical solutions may be easier to obtain in this case than their VA counterparts. 

The knowledge of the exact shape of the solitons may be essential for experiments, as the available size of MM 
samples, as well as the size of PCs, may be quite small, hence the transmission distance in the sample may not be long 
enough to allow self-shaping of the input signal into the soliton. Therefore, to observe and use the soliton propagation, 
it may be necessary to prepare an input pulse (in the case of the MM) or spatial beam (to be coupled into a PC) 
whose shape is as close as possible to that of the stable soliton, thus minimizing detrimental shape oscillations in the 
course of the transmission. 



C. Stability regions 

The stability of all the stationary solitons was tested by means of direct simulations of the propagation within the 
framework of Eqs. ((T|) and adding small perturbations to the initial profile. As a result, it has been concluded 
that all the solitons which could be found in the stationary form, as reported above, are stable, including those with 
extended tails, which were found at small positive and negative values of Z?2- Moreover, it has been concluded that 
the propagation distance necessary for the self-healing of weakly perturbed solitons and their relaxation back to the 
unperturbed shape was always essentially smaller than z = 10 (in the present notation). Examples of the stable 
evolution of solitons with complex shapes, which feature strong tails or the split-peak structure, are displayed below 
in Fig. see also other relevant examples in Figs. l5l and [TOTa). 

Thus, the stability region for the bimodal solitons is supposed to be identical to their existence area. For fixed 
material constants, D2 and c (actually, in this section we consider the case of c = 0), soliton families depend on two 
intrinsic parameters, that may be identified as their wavenumbers K u and K v , or, more conveniently, as the total 
energy, E, and the wavenumber ratio, R = K v /K u . An adequate rendition of the existence areas for the solitons of 
both types that were defined above, i.e., single- and double-hump ones, is provided by plotting the areas in the plane 
of (E,D 2 ) at fixed values of R. In Fig. 01 we display overlapping existence regions of both soliton species for two 
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(a) (b) 

FIG. 4: Borders of the regions of the existence (and, presumably, stability) for the single-peak and double-peak solitons, in 
the model with c = (zero group-velocity mismatch) are shown by dashed-dotted and continuous lines, respectively. Panels 
(a) and (b) correspond to fixed values of the wavenumber ratio: R = K v /K u — 2 in (a), and R — 1.2 in (b). The overlap of 
the stability areas for both species of the solitons demonstrates the bistability of the system. In parts of the existence regions 
located at D2 < the ^-component of the solutions, strictly speaking, cannot take a localized form. In fact, the bottom 
existence border separates the solutions which are "practically localized" and those with conspicuous tails, that do not vanish 
in the entire integration domain. For that reason, extension of the domain may lead to slight expansion of the existence region. 
There is no right border for the existence of the solitons, as they persist to indefinitely large values of the energy. 



characteristic values of the wavenumber ratio, R — 2 and 1.2. The large region of the bistability is evident in this 
figure. 

As one approaches the upper existence border for the solitons of either type, increasing D2 at a fixed value of total 
energy E, the energy ratio, E v /E u (see Eq. (J5J), decreases, vanishing at the upper border. Above this border, only 
single-component solitons exist, with v — (of course, these are regular single-peak NLS solitons). The opposite 
trend is observed with the decrease of D% at fixed E: as one approaches the lower existence border, the it-component 
of the soliton tends to vanish, along with its energy share, while the ^-component develops strong tails, especially in 
parts of the existence region located at D2 < 0. It is easy to explain these trends, looking at the Hamiltonian of the 
present model (with c = 0), 



H = 



' (1/2) (K| 2 + D 2 \v T \ 2 + M 4 + H 4 ) + 2\uv\ 2 dr. (8) 



Indeed, the minimization of the gradient part of the Hamiltonian at fixed E is obviously favored by the attenuation of 
the w-component for large D2, and by the decay of the u-component for small Z?2- Another feature observed in Fig. 
[4] is the shrinkage of the existence regions for the double-hump solitons as the wavenumber ratio approaches R = 1 
(as suggested by the comparison of the panels appertaining to R = 2 and R = 1.2). At R — 1 (equal wavenumbers in 
the two component), the system gives rise to single- hump solitons only. 

A mathematically rigorous approach to the stability analysis for the solitons is based on the computation of the 
corresponding eigenvalues, using equations for perturbation modes linearized around the stationary soliton solutions. 
In particular, in Ref. [26| is was concluded that all multi-hump solitons may be unstable, in this sense, in the case 
of D2 = T On the other hand, our direct simulations of the evolution of double- hump solitons did not reveal any 
tangible instability at D2 = 1 (note that line D 2 = 1 goes across the solitons' existence regions in Fig. weakly 
perturbed solitons remain apparently stable in this case, while a strong perturbation excites intrinsic vibrations in 
the soliton, but does not destroy it, see an example in Fig. O A possible explanation to these findings (thanks to 
which the solitons are classified here as stable also at D2 = 1) is that the respective instability eigenvalues may be 
very small. 



FIG. 5: The perturbed evolution of the soliton for D2 ~ 1, when its energy, E = 5, was suddenly increased by 5%. The 
wavenumber ratio of the unperturbed soliton is R = 2. 



D. The split-step method and "tailed" solitons 



As said above, the application of the pseudospectral simulation algorithm readily corroborates the stability of all 
solitons that can be found in the stationary form, including those which feature conspicuous "tails" (see, e.g., Fig. 
[3]). On the other hand, it is relevant to mention that the popular split-step-fast-Fourier-transform (SSFFT) method 
may produce numerical problems in the case of strong "tails". An example of this is shown in Figs. [6] and [3 which 
demonstrate the difference in the numerical stability of the SSFFT algorithm in the absence of well-pronounced tails, 
and in cases when they are present. 

We stress that the evolution of the solitons was simulated by means of the SSFFT without adding any initial 
perturbations, i.e., Figs. [Jja) and (b) display results generated by an intrinsic instability of the algorithm when it is 
applied to the tailed and double-hump solitons. The SSFFT was realized, in these examples, covering the integration 
domain by 2048 points. 

In fact, the solitons which were used as initial conditions in Figs. [5] and [7| are fully stable. To demonstrate this, in 
Fig. [8] we display simulations of the same pair of solitons as in Fig. [Jj but performed by means of the pseudospectral 
algorithm which was outlined above. 



As said above, dissipative effects may be negligible under experimentally relevant conditions in some types of 
artificial optical media, such as PCs and PC fibers, but they constitute an essential ingredient in the description of 
the light transmission through MMs [Tt], EH • In the simplest approximation, the loss is taken into account by means 
of a positive dissipative coefficient, a, added to Eqs. |T]) and @: 



E. Effects of the loss and its compensation 




(9) 
(10) 



iv z + icv T + (l/2)D 2 u TT + (\v\ 2 + 2|u| 2 ) v 



= —lav. 



Obviously, the action of the loss leads to the decrease of the soliton's energy. As seen in Fig. 0] the gradual 
decrease of the total energy should drive the solitons in the direction of the transition from the double-hump shape 
to the single-hump one, i.e., the loss is expected to make the solitons' shape "more primitive". This expectation is 
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FIG. 6: Results of the application of the split-step fast-Fourier-transform (SSFFT) method to the simulation of the evolution 
of a single-peak soliton without well-pronounced "tails", for Di = 0.3 and E — 2, R = E v /E u — 2. In this case, the SSFFT 
algorithm is relatively stable. 




(a) (b) 

FIG. 7: Unstable operation of the SSFFT algorithm simulating the evolution of a single-peak soliton with conspicuous tails 
(a), and a double-peak soliton (b). In both cases, the relative-GVD coefficient in Eq. ([2]) is D2 ~ 0.1, and the initial solitons 
have E — 2 and R = K u /K u = 2. 

corroborated by a typical example of the evolution of the soliton under the action of the loss displayed in Fig. [9(a): 
under the action of the loss, the soliton gradually decays, transiting from the "tailed" shape to a tailless one. 

As mentioned above, a number of works considered possibilities to periodically compensate the loss in MMs by an 
optical gain [1, [l|| H(| . I n Fig- El we display an example of that, assuming that, after passing each interval Az = 1, 
fields u(t) and v(t) are multiplied by a common factor G = exp (aAz), whose value is selected so as to compensate 
the loss in Eqs. © and (|10p . It is seen from the figure that the periodic compensation may readily maintain the 
soliton's tailed shape. 



(a) (b) 
FIG. 8: Simulations of the same initial solitons as in Fig. [7] but carried out by means of the stable pseudospectral method. 




(a) (b) 

FIG. 9: (a) Evolution of a typical soliton, with an originally "tailed" shape, initial total energy E = 5, and Ro = K v /K u = 1.13, 
under the action of the loss with a — 0.05 and D2 = 0.37 in Eqs. © and (|10[) . (b) The same, but with the compensation of 
the loss by means of the instantaneous linear amplification, applied with period Az = 1. 



III. STABILIZATION OF SOLITONS BY MEANS OF THE GROUP- VELOCITY-MISMATCH 

MANAGEMENT 



In this work, we do not aim to find numerically exact complex (chirped) stationary solutions that may exist in the 
presence of the GVM, c ^ (such solutions are known in the \^ model with the GVM between the fundamental 
and second harmonics, which is related to the existence of "walking" solitons 41 1, and also in the system of 
coupled cubic NLS equations describing the co-propagation of orthogonal linearly-polarized modes i n op tical fibers, 
that includes, in addition to the XPM interaction, the nonlinear coupling via the four- wave mixing [42j|). Actually, 
it would be quite difficult to prepare such a pulse with a specific distribution of the internal chirp for the use in the 
experiment. On the other hand, an experimentally relevant problem is to simulate the evolution of the previously 
found stationary solitons in the framework of Eqs. ([1]) and @ that include the GVM term. To this end, in Fig. [TUTa) 
we display an example of the evolution of the soliton in the presence of relatively weak GVM (the particular value, 
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(a) 



(b) 



FIG. 10: Evolution of solitons in the presence of the group- velocity-mismatch term, with c = —0.144 (a) or c = — 1 (b). Other 
parameters are D2 = 0.37 and E — 5, R = K v /K u — 1.13. 

c = —0.144, used in this figure, corresponds to physically relevant parameters as per the above-mentioned lossless 
Drude model 0]). It is seen that, while generating a permanent perturbation of the soliton, the moderate GVM term 
does not destroy it. On the other hand, Fig. [TOTb) demonstrates that stronger GVM, with c = 1, gives rise to a very 
strong perturbation, under which the soliton cannot maintain its integrity, even in an approximate form. 

Because MMs and PCs are artificially built media, it may be quite natural to consider possibilities to engineer 
superstructures in them that can implement the GVM management, i.e., periodic alternation of layers with large 
GVM coefficients of opposite signs. In this way, it turns out easy to stabilize solitons against the locally strong GVM, 
provided that the management period, L, is relatively small in comparison with the effective wavelength, A„ = 2w/K v . 
A typical example of the stabilization is displayed in Fig. 1111 The scheme used in this figure is based on the following 
management map, with the zero average value of the GVM coefficient: 



which repeats with period L. In the case shown in Fig. 1111 L = 0.2, which amounts to A„/4. 

Finally, it is relevant to mention that, except for the case of D-x = 1 in Eq. $Z§, the system as a whole does not 
support the Galilean invariance, therefore moving solitons may be essentially different from the quiescent ones, that 
were considered above. However, the study of moving localized solutions within the framework of Eqs. ([1]) and |(5J) is 
actually tantamount to considering solitons under the action of the GVM effect. 



The objective of this work was to study solitons in a model based on the system of two XPM-coupled NLS equations, 
which describes the co-propagation of two waves in MMs (metamaterials) , if the loss may be neglected, and in PCs 
(photonic crystals). The same system is relevant for ordinary optical fibers, close to the zero-dispersion point. The 
main peculiarity of the system is the small positive or negative value of the GVD coefficient in one equation, while 
the dispersion is fixed to be anomalous in the other. Unlike previously studied systems of nonlinearly coupled NLS 
equations with equal GVD coefficients, which gives rise only to simple stable single-peak solitons with an arbitrary 
ratio of the energy between the components, the present model generates soliton families which feature complex 
shapes, including those with weakly localized oscillatory tails and a double-peak maximum. Existence regions for 
the single- and double-peak solitons have been identified in the parameter planes, demonstrating a broad range of 
the bistability in the system. Above the upper existence border, the two-component solitons degenerate into obvious 
single-component fundamental NLS solitons. Systematic direct simulations demonstrate that solitons of both types 
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FIG. 11: The stabilization of the same soliton as in Fig. [10] provided by the GVM-management scheme based on map (|lip 
with L — 0.2 — 0.25 • (2tv/K v ). 



are stable throughout their existence regions. Under the action of the loss, the decaying solitons tend to evolve towards 
more "primitive" shapes, losing their tails. On the other hand, the periodic compensation of the loss helps to maintain 
the original shape. The numerical calculations were performed by means of a properly adjusted pseudospectral method 
(while the usual split-step algorithm may be problematic when applied to simulating the evolution of solitons with 
conspicuous tails). The extended model, including the GVM (group- velocity mismatch) was studied too. It has been 
demonstrated that, while a sufficiently strong GVM term tends to destroy the bimodal solitons, they can be readily 
stabilized by way of the GVM-management scheme. 

This work suggests extensions in several directions. First, a study of collisions between moving solitons may be 
a natural addition to the above analysis. Dispersion management |43( may also be considered in the framework of 
the present model. On the other hand, the model can be made two- and three-dimensional by adding transverse 
coordinates. In that case, a challenging issue would be a possibility to stabilize the respective spatiotemporal solitons 
("light bullets") against the collapse by means of a grating created in the transverse direction(s), as well as by means 
of the respectively modified dispersion-management scheme, cf. Ref. fjij . 
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